Load in libraries

library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(DT)
library(ggpubr)
library(cowplot)
source("scripts/utils/utils.R")
modernacolors()
## <environment: R_GlobalEnv>

Load data

# load ELISpot data
temp <- read_excel("processed_data/ELISpot.xlsx")

dp <- temp[!is.na(temp$counts), ] %>%
  filter(`low cell number (<0.15M)`== "no",
         day != "143")%>%
  mutate(counts = counts + 0.001)

d <-
  dp %>%
  filter(interval != "PBS")%>%
  mutate(animal_id = factor(animal_id),
         interval = factor(interval),
         day = factor(day))

Observed data visualization

dp %>%
  filter(interval != "PBS") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
  ggplot(aes(x = interval, y = counts)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "red",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  facet_wrap(tissue ~ coating, scales = "free") +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "antigen specific \n ASC counts/million cells",
    color = "",
    linetype = "",
    title = "Antigen specific spleen and BM results"
  ) +
  theme(
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "bottom"
  )

Model fitting for spleen data

#-------------------------------------------------------------------------------
# Fit the model for spleen data 
#-------------------------------------------------------------------------------
# boxcox function to fit boxcox regression
boxcox_fit <- function(dat) {
  require(MASS)
  if (length(unique(dat$day)) == 1) {
    out <- boxcox(counts ~ interval, data = dat)
    optimal <- out$x[which.max(out$y)]
    bctran <- make.tran("boxcox", optimal)
    lm_fit <- with(bctran,
                   lm(linkfun(counts) ~ interval, data = dat))
  } else{
    out <- boxcox(counts ~ interval * day, data = dat)
    optimal <- out$x[which.max(out$y)]
    bctran <- make.tran("boxcox", optimal)
    lm_fit <- with(bctran,
                   lm(linkfun(counts) ~ interval * day, data = dat))
  }
  return(list(lm_fit = lm_fit, lambda = optimal))
}

# box cox model fit for spleen data
d_tissue <- subset(d, tissue == "spleen") %>%
  filter(!(coating == "NTD" & interval == "0"))

boxcox_fits <-
  d_tissue %>%
  split(.$coating) %>%
  map( ~ boxcox_fit(.x))

lm_fits <-
  purrr::transpose(boxcox_fits)$lm_fit

# obtain optimal lambda for boxcox transformation
lambdas <-
  purrr::transpose(boxcox_fits)$lambda %>% unlist()

assumption_plots <-
  imap(
    lm_fits ,
    ~ assumption_check(.x) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0("Spleen ", .y, "-specific model residual diagnostics "),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )

# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
  emmeans(lm_grid, specs = ~ interval, type = "response") %>%
  as_tibble()

Residual diagnostics for spleen data

ggarrange(plotlist = assumption_plots)

Figure 4A. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells

#-------------------------------------------------------------------------------
# Figure 4A. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
# fold change figure
foldchange1 <-
  emmt_dat %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(facet = "1 week post-boost") %>%
  ggplot(aes(x = interval , y = response, color = interval)) +
  geom_point(position = position_dodge(width = 0.5), size = 2)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 0.5),
    size = 1
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  #ggtitle("S2P specific ASCs") +
  facet_grid( ~ facet) +
  labs(x = "", y = "Estimated counts/million cells") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(6, 7)],
             color = "darkcyan",
             lty = "dashed") +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

# boxplot
box1 <-
  dp %>%
  filter(interval != "PBS", tissue == "spleen", coating == "S2P") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(facet = "1 week post-boost") %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(factor(day), "day 64" = "64")) %>%
  ggplot(aes(x = interval, y = counts, color = interval)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1500)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "grey",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  facet_grid( ~ facet) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "Counts/million cells",
    color = "",
    linetype = "",
    title = "S2P specific ASCs"
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
    )
  ) +
  theme(
    axis.text = element_text(size = 12),
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "None",
  )


cowplot::plot_grid(
  plotlist = list(box1, foldchange1),
  nrow = 2,
  ncol = 1,
  labels = c("A", "B")
)

Supplementary Table 3. Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing intervals.

#-------------------------------------------------------------------------------
# Supplementary Table 3
#-------------------------------------------------------------------------------
emm_elispot <-
  emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
          specs = ~ interval,
          type = "response")

set.seed(100)
sum_elispot <-
  contrast(
    regrid(emm_elispot, transform = "log10"),
    "pairwise",
    infer = TRUE,
    adjust = "mvt",
    ratio = FALSE,
    level = .95
  ) %>% as_tibble() %>%
  bind_rows()

sum_elispot  %>%
  as_tibble() %>%
  dplyr::select(contrast, estimate, p.value) %>%
  filter(p.value < 0.05) %>%
  mutate(estimate = 10 ^ estimate) %>%
  mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value"
  ) %>%
  DT::datatable(caption = "Statistical comparisons of S2-P–specific ASCs 1 week following dose 2 between mRNA-1273 (10 μg) dosing interval") %>%
  formatSignif(columns = c('Fold change', 'Adjusted P-value'),
               digits = 7)
# Supplementary Table 3 visualization
t3v <-
  sum_elispot %>%
  as_tibble() %>%
  separate(contrast, c("a", "b"), sep = " - ") %>%
  mutate(a = str_remove_all(a , "interval"),
         b = str_remove_all(b , "interval")) %>%
  mutate(
    a = case_when(
      str_detect(a, "1") ~ paste(a, "week"),
      str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
      TRUE ~ paste(a, "weeks")
    ),
    b = case_when(
      str_detect(b, "1") ~ paste(b, "week"),
      str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
      TRUE ~ paste(b, "weeks")
    )
  ) %>%
  mutate(contrast = paste(a, "/", b)) %>%
  mutate(day = "1 week post-boost)") %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
  ) %>%
  mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
  ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
  geom_point(size = 3)  +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = .1,
                #position = position_dodge(width = 2),
                size = 1) +
  scale_color_manual(
    values = c(
      "#ABCFE5",
      "#93C4DE",
      "#78B5D8",
      "#60A6D1",
      "#4A97C9",
      "#3787C0",
      "#2575B7",
      "#1664AB",
      "#09539D",
      "#084185"
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
                     labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
  facet_grid( ~ day) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  labs(title = "S2-P-specific ASCs",
       y = "Estimated Fold-Change",
       x = "Prime-boost dosing interval") +
  coord_flip() +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    strip.text = element_text(size = 10),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

Model fitting for BM data

#-------------------------------------------------------------------------------
# Fit the model for BM data
#-------------------------------------------------------------------------------
# BM
d_bm <- subset(d, tissue == "BM")

boxcox_fits <-
  d_bm %>%
  split(.$coating) %>%
  map(~ boxcox_fit(.x))

lm_fits <-
  purrr::transpose(boxcox_fits)$lm_fit

# obtain emmeans
lm_grid <- ref_grid(lm_fits$S2P, cov.keep = c("interval", "day"))
emmt_dat <-
  emmeans(lm_grid, specs = ~ interval | day, type = "response") %>%
  as_tibble()

Residual diagnostics for BM data

assumption_plots <-
  imap(
    lm_fits ,
    ~ assumption_check(.x) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0("BM ", .y, "-specific model residual diagnostics "),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )
ggarrange(plotlist = assumption_plots)

Figure 4B. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells

#-------------------------------------------------------------------------------
# Figure 4B. Spike-specific Antibody Secreting Cells and Long-Lived Plasma Cells.
#-------------------------------------------------------------------------------
foldchange2 <-
  emmt_dat %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(
    factor(day),
    "4 week post-boost" = "87",
    "24 week post-boost" = "226"
  )) %>%
  mutate(day = factor(day, c("4 week post-boost" , "24 week post-boost"))) %>%
  ggplot(aes(x = interval , y = response, color = interval)) +
  geom_point(position = position_dodge(width = 0.5), size = 2)  +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    width = .1,
    position = position_dodge(width = 0.5),
    size = 1
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  labs(x = "", y = "Estimated counts/million cells") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(6)],
             color = "darkcyan",
             lty = "dashed") +
  geom_hline(yintercept = emmt_dat$lower.CL[c(13)],
             color = "darkcyan",
             lty = "dashed") +
  facet_grid(~ day) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme(
    axis.text = element_text(size = 12),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_blank(),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

box2 <-
  dp %>%
  filter(interval != "PBS", tissue == "BM", coating == "S2P") %>%
  mutate(interval = factor(interval, levels = c("PBS", 0:8))) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  mutate(day = forcats::fct_recode(
    factor(day),
    "4 week post-boost" = "87",
    "24 week post-boost" = "226"
  )) %>%
  mutate(day = factor(day,
                      levels = c("4 week post-boost", "24 week post-boost"))) %>%
  ggplot(aes(x = interval, y = counts, color = interval)) +
  geom_boxplot(outlier.shape = NA) +
  scale_shape_manual(values = c(16, 5)) +
  scale_y_log10(breaks = c(100, 500, 1000, 1600)) +
  geom_jitter(aes(x = interval, y = counts)) +
  stat_summary(
    fun.min = mean,
    fun = mean,
    fun.max = mean,
    color = "grey",
    width = .5,
    geom = "errorbar",
    size = 1
  ) +
  scale_color_manual(
    values =   c(
      modernadarkgray,
      modernadarkblue2,
      modernaorange,
      modernapurple,
      modernagreen,
      modernablue,
      modernared
      
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  labs(
    x = "",
    y = "Counts/million cells",
    color = "",
    linetype = "",
    title = "S2-P-specific LLPCs"
  ) +
  facet_grid(~ day) +
  theme(
    axis.text = element_text(size = 12),
    title = element_text(size = 14),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    ),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.position = "None",
    plot.title = element_text(size = 14, face = "bold")
  )

cowplot::plot_grid(
  box1,
  box2,
  foldchange1,
  foldchange2,
  labels = c("A", "C", "B", "D"),
  align = 'hv',
  rel_widths = c(1, 2),
  rel_heights = c(1, 1)
)

#ggsave(filename = "figures_tables/ELISpot/ELISpot_boxplot_F4.png", width = 14, height=10)
#ggsave(filename = "figures_tables/ELISpot/ELISpot_boxplot_F4.pdf", width = 14, height=10)

Supplementary Table 4. Statistical comparisons of S2-P–specific LLPCs 4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273 (10 μg) dosing intervals.

#-------------------------------------------------------------------------------
# Supplementary Table 4
#-------------------------------------------------------------------------------
emm_elispot_BM <-
  emmeans(ref_grid(lm_fits$S2P, cov.keep = c("interval", "day")),
          specs = ~ interval |
            day,
          type = "response")

set.seed(100)
emm_elispot_BM <-
  contrast(
    regrid(emm_elispot_BM , transform = "log10"),
    "pairwise",
    infer = TRUE,
    adjust = "mvt",
    ratio = FALSE,
    level = .95
  ) %>% as_tibble() %>%
  bind_rows()

emm_elispot_BM %>%
  as_tibble() %>%
  dplyr::select(day, contrast, estimate, p.value) %>%
  filter(p.value < 0.05) %>%
  mutate(day = as.numeric(as.character(day))) %>%
  arrange(day) %>%
  mutate(estimate = 10 ^ estimate) %>%
  mutate(contrast = str_replace(contrast, " - ", " / ")) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
    "Day" = "day"
  ) %>%
  DT::datatable(
    caption = "Statistical comparisons of S2-P–specific LLPCs 4 weeks following
                dose 2 (Day 87) and 24 weeks following dose 2 (Day 227) between mRNA-1273
                (10 μg) dosing intervals"
  ) %>%
  formatSignif(columns = c('Fold change', 'Adjusted P-value'),
               digits = 7)

Supplementary Figure 5. Statistical comparisons of S2-P-specific Antibody Secreting Cells and Long-Lived Plasma Cells between mRNA-1273 (10 µg) dosing intervals(Supplementary Table 3 and 4 visualization)

t4v <-
  emm_elispot_BM %>%
  as_tibble() %>%
  separate(contrast, c("a", "b"), sep = " - ") %>%
  mutate(a = str_remove_all(a , "interval"),
         b = str_remove_all(b , "interval")) %>%
  mutate(
    a = case_when(
      str_detect(a, "1") ~ paste(a, "week"),
      str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
      TRUE ~ paste(a, "weeks")
    ),
    b = case_when(
      str_detect(b, "1") ~ paste(b, "week"),
      str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
      TRUE ~ paste(b, "weeks")
    )
  ) %>%
  mutate(contrast = paste(a, "/", b)) %>%
  mutate(day = factor(day, levels = c(87, 226))) %>%
  mutate(day = forcats::fct_recode(
    day,
    "4 week post-boost" = "87",
    "24 week post-boost" = "226",
  )) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "estimate",
    "Adjusted P-value" = "p.value",
  ) %>%
  mutate(`Pairwise comparison` = str_remove_all(`Pairwise comparison` , "interval")) %>%
  ggplot(aes(x = `Pairwise comparison` , y = `Fold change`)) +
  geom_point(size = 3)  +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = .1,
                #position = position_dodge(width = 2),
                size = 1) +
  scale_color_manual(
    values = c(
      "#ABCFE5",
      "#93C4DE",
      "#78B5D8",
      "#60A6D1",
      "#4A97C9",
      "#3787C0",
      "#2575B7",
      "#1664AB",
      "#09539D",
      "#084185"
    )
  ) +
  theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
  scale_y_continuous(breaks = log10(c(0.1, 0.2, 0.5, 1, 2, 3)),
                     labels = c("1/10", "1/5", "1/2", 1, 2, 3)) +
  geom_hline(yintercept = 0,
             color = "darkcyan",
             lty = "dashed") +
  facet_grid(~ day) +
  labs(title = "S2-P-specific LLPC",
       # subtitle = "4 weeks following dose 2 (Day 87) and 24 weeks following dose 2 (Day 227)
       # between mRNA-1273 (10 μg) dosing intervals",
       y = "Estimated Fold-Change",
       x = "Prime-boost dosing interval") +
  coord_flip() +
  theme(
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 12),
    strip.text = element_text(size = 10),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    plot.title = element_text(size = 14, face = "bold"),
    panel.spacing = unit(0, "lines")
  )

cowplot::plot_grid(
  plotlist = list(t3v, t4v),
  labels = c("A", "B"),
  nrow = 1,
  ncol = 2,
  rel_widths = c(0.5, 1)
)

ggsave(filename = "figures_tables/ELISpot/emmeans_BM_spleen.png",
       width = 16,
       height = 8)
ggsave(filename = "figures_tables/ELISpot/emmeans_BM_spleen.pdf",
       width = 16,
       height = 8)

Session information

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MASS_7.3-51.5   cowplot_1.1.1   ggpubr_0.4.0    DT_0.24        
##  [5] lme4_1.1-28     Matrix_1.2-18   ggeffects_1.1.3 emmeans_1.7.5  
##  [9] mgcv_1.8-31     nlme_3.1-144    readxl_1.4.0    forcats_0.5.1  
## [13] stringr_1.4.0   dplyr_1.0.9     purrr_0.3.4     readr_2.1.2    
## [17] tidyr_1.2.0     tibble_3.1.8    ggplot2_3.3.6   tidyverse_1.3.2
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.2            lubridate_1.8.0     httr_1.4.3         
##  [4] tools_3.6.3         backports_1.4.1     bslib_0.4.0        
##  [7] utf8_1.2.2          R6_2.5.1            DBI_1.1.3          
## [10] colorspace_2.0-3    withr_2.5.0         tidyselect_1.1.2   
## [13] compiler_3.6.3      textshaping_0.3.6   cli_3.3.0          
## [16] rvest_1.0.2         xml2_1.3.3          labeling_0.4.2     
## [19] sass_0.4.2          scales_1.2.0        mvtnorm_1.1-3      
## [22] systemfonts_1.0.2   digest_0.6.29       minqa_1.2.4        
## [25] rmarkdown_2.14      pkgconfig_2.0.3     htmltools_0.5.3    
## [28] highr_0.9           dbplyr_2.2.1        fastmap_1.1.0      
## [31] htmlwidgets_1.5.4   rlang_1.0.4         rstudioapi_0.13    
## [34] jquerylib_0.1.4     generics_0.1.3      farver_2.1.1       
## [37] jsonlite_1.8.0      crosstalk_1.2.0     car_3.1-0          
## [40] googlesheets4_1.0.1 magrittr_2.0.3      Rcpp_1.0.9         
## [43] munsell_0.5.0       fansi_1.0.3         abind_1.4-5        
## [46] lifecycle_1.0.1     stringi_1.7.8       yaml_2.3.5         
## [49] carData_3.0-5       grid_3.6.3          crayon_1.5.1       
## [52] lattice_0.20-38     haven_2.5.0         splines_3.6.3      
## [55] hms_1.1.1           knitr_1.39          pillar_1.8.0       
## [58] boot_1.3-24         estimability_1.4.1  ggsignif_0.6.3     
## [61] reprex_2.0.1        glue_1.6.2          evaluate_0.16      
## [64] modelr_0.1.8        vctrs_0.4.1         nloptr_1.2.2.3     
## [67] tzdb_0.3.0          cellranger_1.1.0    gtable_0.3.0       
## [70] assertthat_0.2.1    cachem_1.0.6        xfun_0.32          
## [73] xtable_1.8-4        broom_1.0.0         coda_0.19-4        
## [76] rstatix_0.7.0       ragg_1.1.3          googledrive_2.0.0  
## [79] gargle_1.2.0        ellipsis_0.3.2